System of elastic hard spheres which mimics the transport properties of a granular gas 
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The prototype model of a fluidized granular system is a gas of inelastic hard spheres (IHS) with 
a constant coefficient of normal restitution a. Using a kinetic theory description we investigate the 
two basic ingredients that a model of elastic hard spheres (EHS) must have in order to mimic the 
most relevant transport properties of the underlying IHS gas. First, the EHS gas is assumed to be 
subject to the action of an effective drag force with a friction constant equal to half the cooling 
rate of the IHS gas, the latter being evaluated in the local equilibrium approximation for simplicity. 
Second, the collision rate of the EHS gas is reduced by a factor ~(1 + a), relative to that of the 
IHS gas. Comparison between the respective Navier-Stokes transport coefficients shows that the 
EHS model reproduces almost perfectly the self-diffusion coefficient and reasonably well the two 
transport coefficients defining the heat flux, the shear viscosity being reproduced within a deviation 
less than 14% (for a > 0.5). Moreover, the EHS model is seen to agree with the fundamental 
collision integrals of inelastic mixtures and dense gases. The approximate equivalence between IHS 
and EHS is used to propose kinetic models for inelastic collisions as simple extensions of known 
kinetic models for elastic collisions. 

PACS numbers: 45.70.Mg, 05.20.Dd, 05.60.-k, 51.10.+y 



I. INTRODUCTION 

As is well known, the prototype model of a granular 
fluid under conditions of rapid flow consists of a gas of 
(smooth) inelastic hard spheres (IHS) characterized by 
a constant coefficient of normal restitution a < 1 P, 0- 
When two particles moving with velocities v and vi col- 
lide inelastically, they emerge after collision with veloci- 
ties v' and v[, respectively, given by 



v" and v" giving rise to post-collisional velocities v and 
Vi: 

v"=v ^—( g .a)a, V »= V1 + — — (g-S)£, (1.3) 
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where g = v — vi is now the post-collisional relative 
velocity. 

The inelasticity of collisions contributes to a decrease 
of the granular temperature T{t) (proportional to the 
mean kinetic energy per particle in the Lagrangian 
frame), i.e., 



Here, a is a unit vector directed along the centers of the 
two colliding spheres at contact and g = v — V! is the 
pre-collisional relative velocity. The collision rule (|l.lf) 
conserves momentum but energy is decreased by a factor 
proportional to the degree of inelasticity 1 — a 2 , namely 



dT 
~dt 



(1.4) 
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Stated differently, while the component of the relative 
velocity orthogonal to <x does not change upon collision, 
the magnitude of the parallel component decreases by a 
factor a: g' ■ a = — cvg • <r, where g' = v' — vi is the post- 
collisional relative velocity. From Eq. (|1.1|1 it is straight- 
forward to get the pre-collisional or restituting velocities 



where £ ~ v(\ — a 2 ) is the cooling rate, v being an ef- 
fective collision frequency. Equation l|1.4|) implies that in 
order to reach a steady state an external energy input is 
needed to compensate for the collisional cooling. 

In the case of a gas of elastic hard spheres (EHS), en- 
ergy is conserved by collisions. However, a cooling effect 
can be generated by the application of a drag force pro- 
portional to the particle velocity, i.e., 



dT 
~dt 
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where 7 is the friction coefficient @- Obviously, the 
choice 7 = makes the drag force produce the same 
cooling effect on the EHS system as inelasticity does on 
the IHS one. As a consequence, at a macroscopic level 
of description, the hydrodynamic balance equations of 



mass, momentum, and energy for the IHS gas are (for- 
mally) identical to those for the frictional EHS gas: 



Inelastic hard spheres Elastic hard spheres (+Friction) 
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Here, D t = dt + u ■ V is the material time derivative, d 
is the dimensionality of the system, m is the mass of a 
sphere, n is the number density, u is the flow velocity, 
T is the granular temperature, P is the pressure tensor, 
and q is the heat flux. The right-hand side term in Eq. 
(|1.8|) comes from Eq. I|1.4|l in the case of IHS, whereas it 
comes from Eq. (|1.5|l (with 7 = i£) in the case of EHS. 

Despite the trivial common structure of the macro- 
scopic balance equations i|1.6|) - (|1.8|) for the free IHS and 
the driven EHS gases, the underlying microscopic dy- 
namics is physically quite different in both systems: in 
the IHS gas (i) each colliding pair loses energy upon col- 
lision but (ii) all the particles move freely between two 
successive collisions; in the EHS case, however, (i) energy 
is conserved by collisions but (ii) all the particles lose 
energy between collisions due to the action of the drag 
force. Therefore, during a certain small time step, only a 
small fraction of particles is responsible for the cooling of 
the system in the IHS case, whereas all the particles con- 
tribute to the cooling in the EHS case. These differences 
are sketched in Fig.^ In principle, there is no reason to 
expect that the relevant nonequilibrium physical proper- 
ties (e.g., the one-particle velocity distribution function) 
are similar for IHS and frictional EHS under the same 
conditions. 

Let us consider, for instance, the homogeneous cooling 
state Q . In that case, the energy balance equation (|1.8|) 
becomes 



dT hc (t) 
dt 



= -Chc(t)T hc (t). 



The solution to Eq. lfO|) is Haff 's law: 
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Thc(0) 
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where we have taken into account that Chc(i) oc T h( ( (t). 
The above cooling law is valid for the homogeneous cool- 
ing state of both IHS and frictional EHS, provided that 
in the latter one considers a time-dependent friction co- 
efficient 7 cx T 1 / 2 and 7(0) = Chc(0). The homogeneous 
Boltzmann equation admits in both cases a scaling solu- 
tion of the form 



/ hc (v, t) = n [m/2 T hc (t)f 2 Q c (c(t)) , 
c(t) = v/V^p. 



(1.11) 
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FIG. 1: Sketch of the microscopic dynamics of inelastic hard 
spheres (IHS, left panel) and elastic hard spheres under a 
friction force (EHS, right panel). At t — 0~ both systems 
are prepared in the same microstate and particles A and B 
are about to collide. In the IHS system, particles A and B 
recede immediately after collision (at t — + ) with a relative 
velocity smaller than the pre-collision one, so the mean kinetic 
energy decreases, -E(0 + ) < E(0~). During the time step St 
no collision takes place and the particles move ballistically, 
so E(8t) — E(0 + ). In the EHS system, the collision between 
particles A and B is elastic, so E(0 + ) — E(0~). However, 
during the time step St all the particles feel the drag force 
and, consequently, E(St) < _E(0 + ). Note that at t = St the 
respective microstates in the IHS and EHS systems are dif- 
ferent, even if the energy loss -E(0~) — E(St) is the same. 



On the other hand, while the distribution function is 
a Gaussian /h c (c) = ir~ d / 2 e~ c ~ for EHS deviations 

from a Gaussian are present in the case of ihs ma. 

These deviations are measured, for instance, by a nonzero 
fourth cumulant (or kurtosis), 



(l 2 



d{d + 2) 



(1.12) 
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and by an overpopulated high-energy tail f^ c {c) ~ e 
In the case of a gas heated by a white-noise forcing, the 
steady-state distribution function is again a Gaussian for 
frictional EHS, while a 2 ^ and /*(c) ~ e~ Ac ' i/2 for IHS 

HQ. 

The above two examples are sufficient to illustrate 
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that, obviously, the IHS and driven EHS systems are 
not strictly equivalent. On the other hand, the differ- 
ences between the homogeneous solutions for IHS and 
EHS are not quantitatively important in the domain of 
thermal speeds (for instance, | «2 1 S 0-02 for IHS with 
a > 0.7). Therefore, it is still possible that both systems 
exhibit comparable departures from equilibrium in inho- 
mogeneous states, in which case transport of momentum 
and/or energy are the most relevant phenomena. As a 
matter of fact, one of the most distinctive features of 
granular gases, namely the clustering instability, has a 
hydrodynamic origin 0, , and so it must also appear 
in a gas of EHS under the action of a drag force with 
a state-dependent friction coefficient 7 proportional to a 
characteristic local collision frequency. 

The aim of this paper is to investigate to what extent a 
frictional EHS gas can "disguise" as an IHS gas in what 
concerns the transport properties of mass, momentum, 
and/or energy, which are dominated by the thermal do- 
main of the velocity distribution function. If that were 
the case, the significant body of work already available 
for the kinetic theory of normal gases could be exploited 
to provide a practical tool for granular gases. A prelimi- 
nary report of this work has been given in Ref. [8j- In Sec. 
ITTI we construct a minimal model of EHS in the frame- 
work of the Boltzmann equation that intends to capture 
the basic properties of the Boltzmann equation for IHS 
with a given coefficient of restitution a. First, as dis- 
cussed above, the EHS system is assumed to be under 
the influence of a drag force. By equating the right-hand 
sides of Eqs. (|1.4f) and (|1.5fl . the friction coefficient 7 of 
the EHS gas should be chosen as half the cooling rate 
£(a) of the true IHS gas. However, this is not practical 
since £ is a functional of the nonequilibrium velocity dis- 
tribution of IHS and we want an autonomous Boltzmann 
equation for EHS, i.e., an equation that does not require 
the previous knowledge of the solution of the Boltzmann 
equation for IHS. One could autonomously define 7 as 
the same functional of the EHS distribution as ( is of the 
IHS distribution, but this would also result in a too com- 
plicated model. For these reasons, we choose 7 = i£o(a), 
where Co is the cooling rate in the local equilibrium ap- 
proximation. This is consistent with the fact that the 
velocity distribution function of EHS is a Gaussian in 
the homogeneous cooling state, as well as in the steady 
state driven by a white- noise thermostat. The price to be 
paid by the simple choice 7 = i(o is that the equivalence 
between the respective sources of cooling, Eqs. (II .4f) and 
(|1.5|) . is only approximate. As a second ingredient of the 
model, the collision rate of the EHS gas, relative to that 
of the IHS gas, defines a dimensionless parameter (3 that 
can be freely chosen to optimize the agreement with the 
IHS properties. In order to have a clue on an optimal 
choice of (3, the Navier-Stokes transport coefficients of 
both systems are compared in Sec. IIIII As a compro- 
mise between simplicity and accuracy, we simply choose 
(3{a) — ^(1 + a). The (approximate) mapping EHS^IHS 
allows one to extend directly to granular gases those ki- 



netic models originally proposed for conventional gases 
[3, El- This is illustrated in Sec. II VI for the Bha tnag ar- 
Gross-Krook and the ellipsoidal statistical |l2j | ki- 
netic models. The extension of the approximate equiva- 
lence between inelastic and (frictional) elastic particles to 
mixtures, dense gases, and Maxwell models is discussed 
in Appendices lA"HCl respectively. In particular, the ap- 
plications to mixtures and to dense gases reinforce the 
choice (3(a) = ^(1 + a). The paper ends with some con- 
cluding remarks in Sec. [V] 

II. MODEL OF FRICTIONAL ELASTIC HARD 
SPHERES 

A. Basic properties of the Boltzmann equation for 
inelastic hard spheres 

The Boltzmann equation for a gas of inelastic hard 
spheres (IHS) is EE E3, IHl 

(9 t +v-V)/ = /">[/,/], (2.1) 

where /(r, v;t) is the one-particle velocity distribution 
function and j( a > [/, /] is the Boltzmann collision opera- 
tor 

J {a) [fj] = o- d - 1 J d Vl J d£e( g -ao(g£) 

x [a- 2 /(v")/K)-/(v)/( Vl )], (2.2) 

where the explicit dependence of / on r and t has been 
omitted. In Eq. 12.2fl . a is the diameter of a sphere and 
O is the Heaviside step function. The pre-collisional or 
restituting velocities v" and v" are given by Eq. (|1.3ll . 
Of course, the collision operator for elastic hard spheres 
(EHS), J (1) [/,/], is obtained from Eqs. (EHJ and lfO|) 
by simply setting a = 1. 

The first d + 2 moments of the velocity distribution 
function define the number density 

n(r,t) = J dv/(r,v;i), (2.3) 

the nonequilibrium flow velocity 

u(r,t) = — |-r f dvv/(r,v;t) = (v), (2.4) 
n{T,t) J 

and the granular temperature 

(2.5) 

where V(r, t) = v — u(r, t) is the peculiar velocity. The 
basic properties of [/, /] are those that determine 
the form of the macroscopic balance equations for mass, 
momentum, and energy, 

/ dv mv J {a) [fJ] =( I . (2.6) 



4 



By standard manipulations of the collision operator, the 
cooling rate can be written as [Til Il6| 

cm^i-" 2 ) ^ 1 :^ ^ (2.7) 



4dT(&&) T(r,t) 



where 



(V?a) = ^(77) / dv i / dv 2 |v 1 -v 2 | 3 /(r,v 1 ;t)/(r,v 2 ;t) 

(2.8) 

is the average value of the cube of the relative speed. The 
properties (|2.6|) lead to the balance equations (|1.6|l - (|1.8fl 
with the following kinetic expressions for the pressure 
tensor and the heat flux, 



Pij=m J d-vViVjfty), 
q=| |dvy 2 V/(v). 



(2.9) 



(2.10) 



The cooling rate £ is a nonlinear functional of the dis- 
tribution function / through the average (V^)- As a 
consequence, £ cannot be explicitly evaluated unless the 
Boltzmann equation is solved. Nevertheless, a simple es- 
timate is obtained from Eq. I|2.8|l by replacing the actual 
distribution function / by the local equilibrium distribu- 
tion 

/ (v) = n{m/2-KT) d / 2 exp(-mV 2 /2T). (2.11) 
In that case, 

_ te -^r(^)(!) 3 ' 2 , ,2.2, 

where fl d = 2it d / 2 /T(d/2) is the total solid angle. When 
the approximation 12.12fl is inserted into Eq. (|2.7I) . one 
gets the local equilibrium cooling rate 0, 0] 

Co(r ) t)=Co^o(r,t), Co*-'"" 2 " 2 



Ad 



(1-a'), (2.13) 



where 



d -x(T 
^(d + 2) \m 



1/2 



(2.14) 



is an effective collision frequency. The local equilibrium 
estimate (|2.13|) expresses the cooling rate as a functional 
of / through the local density and temperature only. In 
addition, its dependence on inelasticity is simply £o cx 
1-a 2 . 

It is now convenient to introduce the modified collision 
operator 

[/,/] = J (Q) [/,/]-^-(V/). (2.15) 



By construction, the operator JW [/, /] has the proper- 
ties 



dv j mv ] jW [/,/]= f 
hmV 2 / \ 



as follows from direct evaluation using Eq. H2.6(l . 



(2.16) 



B. The model 

The modified inelastic collision operator J 1 -"' [/, /] 
shares with the elastic collision operator jW [/, /] the 
property of having d + 2 vanishing low velocity mo- 
ments. Of course, [/, /] and [/, /] differ in many 
other aspects. For instance, in the homogeneous cooling 
state the velocity distribution function is the solution to 
•^°^[/hcj /he] = 0| , which differs from a Gaussian, the 
latter being the solution to J^'[f, /] = 0. Moreover, the 
elastic collision operator satisfies the H-theorem, namely 



dv(m/)J«[/,/] <0, 



(2.17) 



while an H-theorem has not been proven for J( Q ) [/, /] . 
Despite these differences, the common properties (I2.16f) 
suggest the possibility that the operators [/, /] and 
[/, /] have a similar behavior in the domain of ther- 
mal speeds (i.e., for V < 2y / 2T/m). This expectation 
can be exploited to propose the approximation 



(2.18) 



where /3(a) is a positive constant to be determined later 
on. Its introduction does not invalidate Eq. I|2.16|l but 
allows us to fine-tune the approximate equivalence be- 
tween both operators. In agreement with the spirit of the 
above discussion we further approximate the true cooling 
rate given by Eq. I|2.7I) by the local equilibrium estimate 
(|2.13|) . In summary, our model consists of the replace- 
ment 

J {a) [fJ] - + |^ • (W) , (2.19) 

so that the Boltzmann equation 12.1f> becomes 



2 dv 



V / = /3J (1) [/,/]. (2.20) 



In this model, the gas of inelastic hard spheres with a 
given coefficient of restitution a is replaced by an "equiv- 
alent" gas of elastic hard spheres subject to the action of 
a drag force Fd rag = —m^V with 7 = ^( . While this 
drag force does not affect the conservation of momentum 
on average, this is not so at a microscopic level. To clar- 
ify this point, let us consider the particles inside a small 
box Sr centered about the point r at time t in a given 
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microstate. Because of the action of the drag force, the 
velocities of those particles change during a short time 
interval St as 



v<(t) -> Vi(t + ft) = Vi(«) - 7(r, t) [Vi(t) - u(r, t)] Ji. 

(2-21) 

Summing over all the <5iV particles inside the box, we 
find that, in general, J2ies r v*(i + <ft) ^ Sie5 r v *O0 sincc 
Sie5r v i(*) 7^ SNu(r,t). On the other hand, the mo- 
mentum is conserved when averaging over all the mi- 

crostates, i.e., (J2ieSr v i( t + St )) = (E Je <5r V *W)> where 
we have taken into acoount that, by definition, u(r,t) = 
(^ ig(5r Vj(i))/(<57V). This averaging process is already 
built in the Boltzmann equation Q2.2U)) . 

The fact that in the friction constant the actual cool- 
ing rate £ of the granular gas is approximated by £q ( or > 
equivalently, (V^) — * (^12)0) i s dictated by simplicity 
since it does not seem necessary to retain the detailed 
functional dependence of Eq. (|2.8(l when on the other 
hand the coarse-grained approximation (|2.18|) is being 
used. In other words, the discrepancies due to ( — > £ 
may be expected to be less important than those associ- 
ated with the approximation (|2.18|1 itself. In any case, if 
one wishes to keep the true cooling rate C i n the model 
(|2.19() , £ must be interpreted as a functional of the solu- 
tion of Eq. H2.2t)|l itself, not as a functional of the solution 
of Eq. H2.ll) . Otherwise, Eq. I|2.20|) would not be an au- 
tonomous equation and it would be necessary to know 
the solution of the IHS Boltzmann equation (|2.1|l before 
dealing with the EHS Boltzmann equation H2.20JI . what is 
not only impractical but artificially complicated as well. 

In Eq. (|2.2|) the coefficient of restitution a appears both 
explicitly (by the factor a~ 2 inside the collision integral) 
and implicitly [through the collision rule l|1.3[l ]. In con- 
trast, in the model H2.19[l a appears only explicitly, as 
well as outside the collision integral, through the approx- 
imate cooling rate (o « 1 - a 2 and the parameter (3(a) 
still to be determined. This simplification can be justi- 
fied as long as one is mainly interested in the gross effects 
of inelasticity on the nonequilibrium velocity distribution 
function, while the fine details might not be captured by 
the model (|2.19|l . However, as seen in the companion pa- 
per 0, the reliability of the EHS model turns out to be 
higher than the expected one. 

In principle, there is no reason to expect that the gas 
of EHS which most efficiently mimics the relevant trans- 
port properties of the granular gas is made of particles 
with a diameter a' equal to the diameter a of the inelas- 
tic spheres. If a' and a were equal, then both systems 
would have the same mean free time but not necessarily 
the same rate of momentum and energy transfer upon 
collisions. The effect associated with a' ^ a is accounted 
for by the parameter (3. Sincc JW[/, /] is defined by set- 
ting a = 1 in Eq. Q2.2)) . it is proportional to cr d_1 , not 
to <j' d l . Therefore, (3J^[f,f] is the collision operator 
of EHS of diameter a' = /jW-iV. Alternatively, Eq. 



(|2.20() can be rewritten as 

dt> + v • V - H ^ • V) / = J« [/, /] , (2.22) 

where t' = (3t and V' = d/dr' with r' = (3r. According 
to Eq. (|2.22|) . the original IHS system is replaced by an 
EHS system with the same diameter a' = a, but with 
a friction constant 7' = Co/2/3 and spatial and temporal 
variables scaled by a factor (3 with respect to those of the 
IHS system. In what follows, we will use for the model 
the form i|2.20[) rather than the form l|2.22l) and will view 
(3 as a correction factor to modify the collision rate of the 
equivalent system of EHS, relative to the collision rate of 
the IHS system. This implies that after a certain common 
time interval At the number of collisions experienced by 
the EHS is in general different from that of the IHS. In the 
next Section we make a definite proposal for (3(a) based 
on the comparison between the EHS and IHS Navier- 
Stokes transport coefficients. 



III. NAVIER— STOKES TRANSPORT 
COEFFICIENTS 

A. Stress tensor and heat flux 

The irreversible momentum and energy transport are 
measured by the stress tensor = Py —pSij (where p = 
nT = d _1 TrP is the hydrostatic pressure) and the heat 
flux q. By an extension of the Chapman-Ens kog method 
fli^ l to the case of inelastic collisions 0, EI H2, , 
one gets the Navier-Stokes constitutive equations 



n, 



77 ( Vjitf + VjUi - -V • uSij 



-AVT - /iVn, 



(3.1) 



(3.2) 



where r\ is the shear viscosity, A is the thermal conductiv- 
ity, and \i is a transport coefficient with no counterpart in 
the elastic case. For IHS, the explicit expressions for the 
transport coefficients in the first Sonine approximation 
are given by 0, 0] 



nT 



1 



V 



nT d + 2l + 2a^ c 

\ — z 

mv 2 z/f - 2C,* ' 



T 2 d + 2 C* + a\ c v* x 
mv 2 K-|C*)K-2C*)' 



(3.3) 



(3.4) 



(3.5) 



In these equations, the effective collision frequency vq is 
defined bv Eq. (|2~Tll . 



4 C 



16(l-a)(l-2a 2 ) 



9 + 24d - a(41 - 8d) + 30a 2 (l - a) 



(3.6) 
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is an estimate of the kurtosis of the velocity distribution 
function [cf. Eq. I|1.12|l ] in the homogeneous cooling state 
HE3, and 



(3.7) 



C*=Co*(l+i6 a " 



is the (reduced) cooling rate in the same state. Moreover, 
, _ /dvD:£( Q )/ D 



vq J dv/ D : D 
3 / 2 



Ad 



l-a + -d (1 + a) 1- 



32 



(3.8) 



is the (reduced) collision frequency associated with the 
shear viscosity, where 



d 



and 



A,(V) ee m ViVj - —6a 



u J dv / S • S 

1 + a 



(3.9) 



d — 1 3 . , „. , . 
— + I g(d + 8)(l-a) 



4 + 5d-3(4-d)a hc 

a 2 



512 



(3.10) 



is the (reduced) collision frequency associated with the 
thermal conductivity, where 



S(V) = 



-V 2 - 



d + 2 



T V. 



(3.11) 



In the first equalities of Eqs. (EH and (|3~TU|) . rep- 
resents the linearization of the collision operator 
around the homogeneous cooling state: 



(3.12) 



In the model (|2.19|l the transport coefficients are for- 
mally given by Eqs. H3.3fl - H3.5fl . except that 



-,hc 





(since the homogeneous cooling state solution is now the 
local equilibrium distribution, i.e., f\ lc = /q) and v * and 
v\ are given by the first equalities of Eqs. (|3.8() and (|3.10|) 
with the replacement 

£ (Q) ^/?£ (1) -y^-V, (3.13) 

where the operator Cy> is defined by Eq. (|3.12|l with 
a = 1 and /i lc = /q. Taking into account the properties 

V • d v D tj = 2Ai, V ■ d v S % = 3Si + {d + 2)TV h (3.14) 



one easily gets 



(3 + C 



o • 



(3.15) 



d - 1 a 3 ,* 



(3.16) 



In summary, the transport coefficients of the EHS model 
in the first Sonine approximation are 



nT 



nT d+ 2 



mu 2 



M = 



T 2 d(d + 2) 



mv Q 2(^-1)^(^/3- |Co)' 



(3.17) 



(3.18) 



(3.19) 



So far, the choice of the parameter f3 remains open. 
To reproduce the main trends in the a-dependence of 
the shear viscosity for IHS, let us equate Eq. I|3.17[) to 
Eq. H3.3fl (with a^ c = for consistency). This yields 



= 



1 + a 



1 



2d 



(1-a) 



/3„. 



(3.20) 



Analogously, the model captures the behavior of the ther- 
mal conductivity if p is obtained by equating Eq. (|3.18|l 
to Eq. (J23J) with a| c = 0. The result is 



(3- 



1 + a 



1 34-^ 



(3x. 



(3.21) 



This expression also optimizes the agreement between 
Eqs. and flH^ . i.e., ^ = j3 X - 



B. Self-diffusion 

If in the homogeneous cooling state a group of tagged 
particles (here labeled with the subscript t) have initially 
a nonuniform density nt, the associated velocity distribu- 
tion function / t obeys the Boltzmann-Lorentz equation 



(d t + v • V) f t = [/ t , / hc ] = (3.22) 

As a consequence, a current j t = Jdvv/ t (v) of tagged 
particles appears opposing the concentration gradient. 
Conservation of the number of tagged particles implies 
the continuity equation 

%it + V ■ j t = 0. (3.23) 

In the limit of weak gradients, one has 

j t = -DVn t , (3.24) 
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where D is the self-diffusion coefficient. By standard ap- 
plication of the Chapman-Enskog method in the first So- 
nine approximation, the self-diffusion coefficient of IHS 
can be derived. The result is [25j 



D = 



T 



1 



1 A* 

2^ 



(3.25) 



where 



D 



v J dv v 2 f 



id 



(3.26) 



is the (reduced) collision frequency associated with the 
self-diffusion coefficient. 

In our EHS model, the Boltzmann-Lorentz equation 
(|3.22l) is replaced by 



d t + v • V 



Cod_ 
2 <9v 



(v - u t ) 



ft = /?J (1) [/t,/o] 



= -/?£«/ t (3.27) 



where we have taken into account that the peculiar veloc- 
ity V appearing on the left-hand side of Eq. I|2.19[) must 
now be understood as v — u t , where u t = jt/«t is the 
mean velocity of the tagged particles, in order to make 
Eq. (|3.27|l consistent with the continuity equation l|3.23|l . 
Similarly to Eq. (|3.13|l . the model implies 



C 



(a) 
BL 



Rr (l) _ Co d 
2 dv 



so that 



fi- 



.d + 2 
2cT 



Q 

2 



(3.28) 



(3.29) 



In addition, the presence of the term proportional to u t 
in Eq. (|3.27|) generates an extra term in the Chapman- 
Enskog method such that £*/2 in the denominator of Eq. 
(|3.25l) is replaced by Q. In summary, the self-diffusion 
coefficient corresponding to the EHS model is (in the first 
Sonine approximation) 



D 



T 



d+2 
2d 



(3.30) 



Identifying Eq. (|3.30(l with Eq. H3.25j) (by setting a 
in the latter), one gets 



fi- 



fin- 



(3.31) 



C. Comparison between the IHS and EHS 
transport coefficients 

We have just seen that the suitable choice for the pa- 
rameter (3(a) under the criterion of optimizing the agree- 
ment with the transport coefficients of IHS is not unique. 




FIG. 2: (Color online) Plot of f3 v ( ), fix (•••)> and f3 D (- 

• - • -) as functions of the coefficient of restitution a in the 
three-dimensional case. 



Depending on the transport property of interest, it may 
be more convenient fi = fi n , fi = fix, fi = (3d, or even a 
different choice. Focusing on these three possibilities, we 
note that (3 V < (3d < fix for d — 2 and d = 3. Moreover, 
fix < 1 for all a in the case d = 3, as shown in Fig. [21 This 
indicates that the diameter of the equivalent EHS system 
is smaller than that of the actual IHS system (a' < a), 
this effect being more pronounced with fi = fi^ than with 
fi = fix- " 

Figure compares the four transport coefficients of 
IHS (d = 3) [Eqs. ipO |l -ipr5 )l and with those of 

the "equivalent" system of EHS [Eqs. (jo"T7|) - lprT3j) and 
IPOU))] wit h the choices fi = fi v [Eq. (3125], fi = fix 
[Eq. (323J], fi = (3d [Eq. (333) ], and fi = 1. It can be 
observed that the EHS shear viscosity combined with the 
choice fi — fi v reproduces almost perfectly the IHS shear 
viscosity. A similar situation occurs in the case of the 
self-diffusion coefficient with the choice fi = (3d- This is 
because the influence of a^ c ^ in the cooling rate (* 
[Eq. (|3.7|l ] and in the collision frequencies v* [Eq. (13. 81 ] 
and v* D [Eq. I|3.26|l ] is very small. On the other hand, 
the deviations of the EHS coefficients A and /i with fi — 
fix from the respective IHS coefficients are much more 
important, essentially due to the explicit dependence on 
a,2 C of the numerators on the right-hand sides of Eqs. 
(|3.4|) and (|3.5|l . since the influence of a^ on is again 
rather weak. 

A natural question that arises is whether there exists 
a common choice for (3(a) that reproduces reasonably 
well the a-dependence of the four transport coefficients. 
Figure |31 shows that fi = (3 V , which is excellent in the 
case of rj, strongly overestimates A, (i, and D. Moreover, 
fi — fix is not clearly superior to (3d in the cases of A and 
/x, whereas it is poorer in reproducing 77 and D. Obvi- 
ously, the naive choice (3 — 1 does not capture well the 
a-dependence of the transport coefficients, especially in 
the case of the shear viscosity. Therefore, we propose 
to take fi — (3d = \{1 + a). With this choice, the fric- 
tional gas of EHS has practically the same self-diffusion 
coefficient as the true IHS gas and similar values for the 
coefficients defining the heat flux. Although the choice 
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It is worth mentioning that the transport coefficients 
obtained from frictional EHS with 7 = 5C0 an d P = 
^(1 + a) exhibit a much better agreement with the IHS 
coefficients than the ones corresponding to the inelastic 
Maxwell model 23] . The relationship between the inelas- 
tic Maxwell model and the (frictional) elastic Maxwell 
model is discussed in Appendix [U] 

Before closing this Section, a comment is in order. The 
transport coefficients 77, A, /i, and D have been consid- 
ered here in the first Sonine approximation, both for IHS 
and EHS, what has allowed us to work with explicit ex- 
pressions. Without prejudicing the degree of reliability of 
the first Sonine approximation, it can be understood as 
a useful tool to probe the structure of the linearized col- 
lision operator through some relevant inner products [see 
the first equalities of Eqs. JOJ, (|3~TU|) . and On 
the other hand, a comparison of the first Sonine approx- 
imation for the transport coefficients with direct simula- 
tion data shows a good agreement for 77 |26l |27| and D 
[2H |2S| , but important discrepancies for high inelasticity 
are present in the cases of A and \x j2fj [2|j . 

IV. KINETIC MODELING 
A. BGK model 



FIG. 3: (Color online) Plot of the (reduced) shear vis- 
cosity ri/(nT/vo), thermal conductivity X/(nT /ravo), trans- 
port coefficient /i/ '(T 2 /mvo), and self-diffusion coefficient 
D/(T/mi/o) for three-dimensional IHS ( — ) and the "equiv- 
alent" system of EHS with = V ( ), P = P\ {•■■), 

P = Pd (- • - • -), and (3 = 1 (-••-•• -) as functions of the 
coefficient of restitution a. Note that in the top panel the IHS 
curve and the EHS curve with P — f3 v are practically indistin- 
guishable. The same happens in the bottom panel between 
the IHS curve and the EHS curve with /3 = Pd- 



/3 = 4(1 + a) underestimates the shear viscosity, this is 
not a serious drawback since this is the transport coeffi- 
cient least sensitive to inelasticity. For instance, the ratio 
between the value of a transport coefficient at a = 0.5 
and the value at a = 1 is about 1.3 in the case of 77, 
while it is about 1.9 and 1.8 in the cases of A and D, 
respectively. At that rather high inelasticity (a = 0.5), 
the EHS transport coefficients with P — i(l + a) differ 
from the IHS ones by 14% (77), 2% (A), 3% (/i), and 0.5% 
(D). 

There are two additional reasons to favor Eq. I|3.31|l . 
First, it is much simpler than Eqs. (|3.20() and (|3.21() and 
does not depend on the dimensionality d. The second 
reason is more compelling. The extension of the model 
(|2.19l) to dilute mixtures and to dense gases is carried out 
in Appendices lAl and FBI respectively. In both cases the 
equation for the collisional transfer of energy imposes the 
choice (I3.31f) as the most natural one, without having to 
resort to the evaluation of transport coefficients. 



The (approximate) mapping IHS<->EHS allows one to 
take advantage of the existence of simple models for EHS 
to extend them straightforwardly to IHS. For instance, 
consider the well-known Bhatnagar-Gross-Krook (BGK) 
model for elastic hard spheres [ill ]: 

J {1) [fJ]^-Mf-fo), (4.1) 

where fo is the local equilibrium distribution (|2.11|) and 
the effective collision frequency vq is usually identified 
with Eq. H2.14f> in order to make the shear viscosity agree 
with that of the Boltzmann equation. Thus, in consis- 
tency with the approximation (|2.19|l . the extension of the 
BGK model to inelastic hard spheres would simply be 

J (a) [f, /] - -PMf ~fo) + ~' (V/) (4.2) 

2 G»V 

with £0 and /3 given by Eqs. I|2.13|l and l|3.31|l . respec- 
tively. In fact, the kinetic model 1)4. 2|) can be seen as a 
simplification of the one already proposed in Ref . ^(| : 

J {a) [f, /] - -PMf - /he) + ^ ■ (V/) , (4.3) 

where here /he represents the local form of the homo- 
geneous cooling state, C is given by Eq. I|3.7|l . and P is 
assumed to be given by Eq. (|3.20|) . 

The linearized collision operator corresponding to the 
BGK model l|4.2|l is simply 

^ _ I iL.V. (4.4) 
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As a consequence, the relation (|3.15(l holds again but 
one has !/£—►/?+ §(q instead of l|3.16f) . Therefore, the 
shear viscosity is given by Eq. H3.17|) . while the thermal 
conductivity and the coefficient jj, are 



nT d + 2 



2 (3-lQ> 



V = 



T 2 d+2 



Co* 



(4.5) 



(4.6) 



which differ from Eqs. I|3.18|l and H3.19[l . respectively, by 
a factor {d — l)/d in front of (3. 

In the case of self-diffusion in the homogeneous cooling 
state, the BGK kinetic equation for tagged particles is 



9 t + v-V-^^--(v-u t ) 
2 av 



ft = -Pvo (ft - ^/o) , 



so the Boltzmann-Lorentz operator becomes 



(4.7) 



4l 5 - P"o 



2 dv 



(4.8) 



and one has v* D — > (3 + §Co instead of (|3.29[l . Thus, 



D = 



T 



mv (3 - ±Co 



(4.9) 



which differs from Eq. 

The inability of the BGK model to reproduce simul- 
taneously the different transport coefficients is already 
present in the elastic case and is the price to be paid by 
the inclusion of a single collision frequency vq. In par- 
ticular, the BGK model yields the value Pr = 1 for the 
Prandtl number Pr = (d + 2)rj/2mX in the elastic limit, 
while the correct Boltzmann value is (in the first Sonine 
approximation) Pr = (d — l)/d. 



B. Ellipsoidal statistical model 

To avoid in part the above limitation, the so-called el- 
lipsoidal statistical (ES) model was proposed about forty 
years ago for elastic particles [3, H2 ■ This kinetic model 
reads 



^ (1) [/,/]^-^oPr(/-/fi), 



(4.10) 



where 



Mv) = n O ^ (det R)-V exp (- ^ R- : w) 

(4.11) 

is an anisotropic Gaussian distribution with the tensor R 
given by 



R = p; ¥ - (i - P] 



(4.12) 



P being the pressure tensor. The ES choice of fn is 
based on information theory arguments. Note that the 
ES model reduces to the conventional BGK model in the 
special case Pr = 1. It is then convenient to consider Pr 
as a free parameter of the model so that the BGK model 
is recovered by formally setting Pr = 1. The reference 
distribution fn has a finite norm provided that R is a 
positive definite matrix, i.e., its eigenvalues must be 
non-negative. From Eq. I|4.12|l . — Pr -1 [p— (1 — Pr)p;], 
where pi are the eigenvalues of the pressure tensor P. 
Since 53i=i Pi = dp, then pi < dp and, consequently, the 
positiveness of implies that Pr > (d— l)/d. The lower 
bound coincides with the physical value of the Prandtl 
number. The first few moments of fn are 



j dv{l,V,mVV}^(v) = {n,0,R}. 



(4.13) 



While in the BGK equation (|4.1|) the reference function 
fo (namely, the local equilibrium distribution) is a func- 
tional of / through its hydrodynamic fields n, u, and T, 
in the ES model fn depends also on the irreversible part 
of the momentum flux. 

When (|4.10() is inserted into (|2.19(l we get our extension 
of the ES kinetic model for IHS: 



J (a) [/,/]^-/3^oPr(/- 



/h) + |^"(V/). (4.14) 

We note that the ES model 1)4. 14|) is different from the 
more detailed Gaussian kinetic model recently proposed 
by Dufty et al. 30), the most important difference being 
that the collision frequency becomes a function of the 
peculiar velocity in the latter Gaussian model. The ki- 
netic model (|4.14|) should also be distinguished from the 
ansatz of an anisotropic Gaussian (or maximum-entropy) 
velocity distribution function introduced by Jenkins and 
Richman |3l| as a means to obtain a closed set of equa- 
tions for the pressure tensor. 

The linearized collision operator associated with (|4.14|) 
is given by (|3.13l) , where now the action of the linearized 
operator in the elastic case is [l(| 



Pr0(v) + ^/o(V)D(V) 

dv' D (V) 0(v') . (4.15) 



This implies that the shear viscosity is given by Eq. 
(|3.17(l . while the thermal conductivity and the fi coef- 
ficient are 



A = 



nT d+2 



1 



2 Pr/3-iCo*' 



T 2 d + 2 



Co 



mv 2 Pr/3(Pr/3- U*)' 



(4.16) 



(4.17) 



respectively. Equations l|4.16|l and (|4.17|l coincide with 
Eqs. ijXTHj) and if one sets Pr = (d - l)/d. Of 
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course, the BGK results (|4.5|l and (|4.6(l are recovered by 
formally setting Pr = 1. 

The natural version of the ES model to the self- d t P xy + (f3v + Co) Pxy + aP yy = 0, (4.24) 



diffusion problem is 

ft = -PVqPt (/ t - ^/ ) • d t Pyy + {pw + Co) Pyy - faop = 0. (4.25) 



9*+v-V--— (v-u t ) 

(4.18) 

Therefore, 



Their steady-state solution is 



D = ^™3i7F> ^ T -^m~QWTor^ (4 - 26) 



mv a Pr/3 - ±Co 

which differs from Eq. 13.30(1 . unless one allows the 

Prandtl number to take the artificial value Pr = (d + yy = _ (4 07) 

2)/2d. nT /3 + Co' 



C. Solution of the BGK and ES models for 
uniform shear flow 

As an illustration of the BGK and ES models extended 
to 1HS, let us analyze their solutions in the case of one 
of the paradigmatic nonequilibrium states, namely the 
uniform (or simple) shear flow. In most of this Subsection 
we will consider the nonlinear ES model l|4.14|l with an 
arbitrary value for Pr, so that Pr = 1 corresponds to the 
BGK model and Pr = (d — l)/d corresponds to the true 
ES model. 

In the uniform shear flow, the density is constant, the 
granular temperature is uniform, and the flow velocity 
has a linear profile u = aj/x, a being the constant shear 
rate. At a more fundamental level, the velocity distri- 
bution function becomes uniform when the velocities are 
referred to the co-moving Lagrangian frame: 

/(r,v;i) = /(V,*), V = v-ayx. (4.20) 

In that case, the ES model kinetic equation reads 

'dvj 2 dV 



dtf~aV y — f~^ — -(Vf)=-f3v Pr(f-f R ). (4.21) 



Multiplying both sides by mViVj and integrating over 
velocity, we get 

d t P i3 + a {8 ix P V j + S ]X P iy ) + CoPij = -(3v (P^ - pSij) , 

(4.22) 

where on the right-hand side we have made use of Eqs. 
(|4.12|) and l|4.13|l . The set of equations l|4.22|l is common 
to the BGK and the ES models since the constant Pr does 
not appear. The structure of Eq. (14.22(1 is also obtained 
from other BGK-like models |32ll33| |. as well as from the 
Boltzmann equation in Grad's approximation |34ll35l |. In 
this latter case, however, the flexibility of accommodating 
the coefficient (3 disappears since by construction it is 
constrained to /3 = f3 n . 

The three independent equations stemming from Eq. 
(ET32|) are 



dtp + CoP ■ 



2a 

d xv 



(4.23) 



p 

r xy 

nT 



2/3 + Co* 



(4.28) 



where in Eq. KM T is an arbitrary reference tempera- 
ture and vq(To) is its associated collision frequency. Ex- 
cept P X y = Pyx, the remaining off-diagonal elements of 
the pressure tensor vanish. In addition, P yy — P zz — 



■ ■ ■ = Pdd, so that P xx = dp 
((4.26(1 can be rewritten as 



(d — I) P yy . Equation 



(4.29) 



which gives the shear rate in units of the steady-state col- 
lision frequency. The rheology of the uniform shear flow 
can be conveniently characterized by a (dimensionless) 
non- Newtonian viscosity coefficient 



Pxy Vq{T) 

nT a 







03 + Co) 5 



(4.30) 



where use has been made of Eqs. ((4.28(1 and ((4.29|) . 

Figure 0| shows the nonlinear shear viscosity ((4.30|) as 
a function of the coefficient of restitution a and of the re- 
duced shear rate a/v$(T) for d = 3 and two choices of (3: 
(3 = P v and /3 = (3d = \{l + a). Numerical data obtained 
from Monte Carlo simulations of the Boltzmann equation 
for IHS H, are also shown. We observe that, except 
perhaps in the quasi-elastic limit, the choice (3 = (3d ex- 
hibits a better global agreement than the choice (3 = (3 n , 
even though the latter is tuned to reproduce the Newto- 
nian shear viscosity (see Fig. Since, as said above, 
Eqs. 1(4.22(1 - 1(4.30(1 with (3 — I3 V are derived from the orig- 
inal Boltzmann equation for IHS in Grad's approxima- 
tion, we remark that the BGK and ES kinetic models 
with the choice (3 = h(l + a) are more accurate than 
Grad's approximation of the Boltzmann equation for this 
particular state. This paradoxical result is partly due 
to the inherently non-Newtonian character of the steady 
uniform shear flow |35] . It is interesting to note that 
for large inelasticity the kinetic models tend to underes- 
timate the shear thinning effect of rj* (see left panel of 
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FIG. 4: (Color online) Steady-state non-Newtonian shear vis- 
cosity for three-dimensional inelastic hard spheres under uni- 
form shear flow as a function of the coefficient of restitution a 
(left panel) and of the reduced shear rate a/vo (right panel). 
The curves are the common predictions of the BGK and ES 
kinetic models, Eqs. (14.2911 and H4.30L with j3 = f3 v (dashed 
lines) and f3 = /3d = (1 + a)/2 (solid lines). The circles are 
simulation results [Isf. 



1.5 nr 




a 



FIG. 5: (Color online) Steady-state fourth- and sixth-degree 
moments, relative to their local equilibrium values, for three- 
dimensional inelastic hard spheres under uniform shear flow 
as functions of the coefficient of restitution a. The curves 
are the predictions of the BGK and ES kinetic models with 
P — Pn (dashed lines) and f3 = /3d = (1 + ot)/2 (solid lines). 
The circles are simulation results pj|. 



Fig. 0} . Since they tend to underestimate the value of 
a/vo(T) as well (not shown), it turns out that the agree- 
ment in the plot rj* versus a/i>o(T) is fairly good even for 
large values of the reduced shear rate (see right panel of 
Fig. © . 

A practical advantage of kinetic models is the pos- 
sibility of obtaining explicitly the velocity distribution 
function. The stationary solution to the kinetic equation 
(14.21(1 can be expressed as 



/(V) = ^PtA-VkCV) 

= /^ Pr / ds exp (-As) f R (V), (4.31) 
Jo 



where the operator A is 



Co- 



A = ^ Pr--Co-<^— -fV--. (4.32) 

The operators V y d/dV x and V ■ d/dV commute. There- 
fore, 

/(V) = /3z/ Pr / ds exp 
Jo 

xf R (e^ s / 2 (V + asV y x)) , (4.33) 
where we have taken into account the properties 

(4.34) 



- (/3*/ Pr--Co 



( e w 2 v) 



(4.35) 



exp I asVy— ) 0(V) = (/> (V + asV v Sc) . 



Since the pressure tensor, and hence the reference dis- 
tribution function f R , are entirely known, Eq. (|4.33|) 
gives the explicit solution. From it one can compute any 
desired velocity moment. Some of those moments are 
derived in Appendix IdI Figure shows (V 4 ) / (V 4 )q and 
(V 6 )/{V 6 ) , where (V 2k ) = (2T/m) k T(k + d/2)/T(d/2) 
is the local equilibrium value, as given by the BGK 
(Pr = 1) and ES (Pr = (d-l)/d)) kinetic models, as 
well as by Monte Carlo simulations of the Boltzmann 
equation for IHS [l^. Even though the ES model is 
more sophisticated than the BGK model, it gives much 
poorer predictions for the fourth- and sixth-degree mo- 
ments than the BGK model. This situation is similar 
to that found in the elastic case j3|J. While this ap- 
pears as a paradoxical result, one must bear in mind 
that the real advantage of the ES model over the BGK 
model occurs when the Prandtl number plays a role, i.e., 
in states where momentum and energy transfer coexist. 
Since in the steady uniform shear flow no energy trans- 
port is present, there is no reason a priori to expect the 
ES model to perform better than the BGK model. On the 
other hand, in the steady Couette flow (where a quasi- 
parabolic temperature profile coexists with a quasi- linear 
velocity profile), both models differ already at the level 
of the non-Newtonian transport properties, the ES model 
presenting in general a better agreement with simulation 
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results 0, 113 • In what concerns the influence of /3, we 
observe in Fig. that (3 = (3 D = (1 + a)/2 is better for 
small or moderate inelasticity, while (3 = (3^ tends to be 
better for large inelasticity, an effect that contrasts with 
that of Fig. H 



V. CONCLUDING REMARKS 

When one aims at understanding the basic properties 
of fluidized granular media, the prototypical model con- 
sists of a system of inelastic hard spheres (IHS) with a 
constant coefficient of restitution a. At a microscopic 
level of description, kinetic theory has proven to be 
an efficient tool to unveil some of the most intriguing 
features of IHS, such as the lack of equilibrium Gibb- 
sian states, high energy overpopulations, clustering in- 
stabilities, Maxwell-demon effects, breakdown of energy 
equipartition, heat fluxes induced by density gradients, 
segregation in mixtures, phase transitions in velocity 
space, . . . The kinetic theory description is based on the 
Boltzmann and Enskog equations for dilute and dense 
granular gases, respectively. Both of them assume the 
molecular chaos hypothesis, according to which the veloc- 
ities of two particles that are about to collide are uncorre- 
cted. Although the degree of validity of this hypothesis 
is much more restricted in inelastic collisions than in the 
case of elastic collisions |38|, the Boltzmann and Enskog 
equations provide useful insights into the peculiar behav- 
ior of granular fluids. 

Mathematically speaking, the Boltzmann equation for 
IHS is much more involved than for elastic hard spheres 
(EHS). As Eqs. (JOJ) and show, the coefficient of 

restitution a appears inside the velocity integral of the 
Boltzmann collision operator in a two-fold way. First, it 
appears explicitly as a factor a~ 2 in front of the gain term 
as a consequence of the properties dv"dv'/ = a _1 dvdvi 
and g" <7 = — a^g-fT. Second, the coefficient of restitu- 
tion appears through the pre-collision velocities v" and 
v" given by the collision rule l|1.3|) . The primary conse- 
quence of a < 1 is the collisional loss of energy, Eqs. (|1.4|) 
and l|2.6l) . which takes place at a cooling rate £ given by 
Eqs. (E3 and (|2~%|) . 

The point we have addressed in this paper is the pro- 
posal of a model of granular gases based on the Boltz- 
mann equation for elastic hard spheres (EHS). In this 
way, the complexities of inelastic collisions for granular 
kinetic theory are represented by a simple two-fold mod- 
ification of the more familiar kinetic theory for elastic 
collisions. Since elastic collisions conserve energy, the 
fundamental ingredient of the model is the introduction 
of some kind of external driving that produces a macro- 
scopic cooling effect similar to the one due to inelastic- 
ity. While the choice of that driving is not unique, the 
most intuitive possibility consists of assuming that the 
EHS gas is under the influence of a (fictitious) drag force 
Fdrag = —7717V. In order to produce the same cooling 
effect than in the true IHS gas, the friction constant must 



be adjusted to be 7 = |C- However, since in general £ 
is a complicated nonlinear functional of the one-particle 
velocity distribution function and we want to keep the 
frictional EHS model as simple as possible, it is prefer- 
able to take 7 = |^0) where Co oc rtT 1//2 (l — a 2 ) is the 
cooling rate in the local equilibrium approximation. The 
second ingredient of the model is subtler. The appli- 
cation of the drag force is not enough for the EHS gas 
to mimic, at least quantitatively, some other collisional 
properties of the IHS gas, apart from the cooling effect. 
In our model we have assumed that the collisions in the 
EHS gas are slowed down by a factor (3 < 1 with respect 
to the IHS gas. In the dilute limit, this can be interpreted 
as assuming that the diameter a' of EHS is smaller than 
the diameter a of IHS, namely a' /a = /J 1 /^ -1 ). Alterna- 
tively, one can interpret /3 as a scaling factor in space and 
time. Comparison between the Navier-Stokes transport 
coefficients derived from the Boltzmann equation (in the 
first Sonine approximation) both for IHS and frictional 
EHS suggests (3 — |(1 + a). Moreover, this choice is 
further supported by the mapping IHS^EHS applied to 
dilute gas mixtures (see Appendix as well as to dense 
single gases (see Appendix |Bj . In the latter case, due 
to the physical separation between two colliding parti- 
cles and the presence of the pair correlation function in 
the Enskog collision operator, the interpretation of (3 in 
terms of different sizes a' =/= a or as a spatio-temporal 
scaling factor is not literally correct. In that case, one 
must understand (3 just as a correction factor modifying 
the collision rate of EHS relative to that of IHS. 

It is obvious that an IHS gas and a gas made of fric- 
tional EHS differ in many respects. For instance, the 
latter admits Maxwellians as uniform solutions (with a 
time-dependent temperature in the case of the homoge- 
neous cooling state or as a stationary one in the case of a 
white-noise forcing) , while the IHS gas typically exhibits 
overpopulated high energy tails. However, even though 
those discrepancies might be qualitatively relevant, they 
are not quantitatively important in the domain of ther- 
mal speeds. In particular, all but the two first distinctive 
features of granular gases listed in the first paragraph of 
this Section can be expected to appear in an EHS gas 
with a state-dependent friction constant 7 cx nT 1 / 2 . 

When the gas is in anisotropic and/or inhomogeneous 
states because of either a transient stage or an instabil- 
ity or the application of boundary conditions, transfers 
of momentum and/or energy are present, so the velocity 
distribution function can deviate strongly from a (local) 
Maxwellian, even for small and moderate velocities. It is 
in those situations where the frictional EHS gas may be 
expected to describe the main transport properties of the 
IHS gas. We have tested this expectation by performing 
Monte Carlo simulations of the Boltzmann equation un- 
der uniform shear flow both for IHS and EHS [|Q1. Th e 
results show that for moderate inelasticity (say a = 0.9) 
one can hardly distinguish the curves representing the 
transient profiles and the steady-state velocity distribu- 
tion functions corresponding to both systems. For larger 



13 



inelasticity (say a = 0.5) the EHS model still captures 
almost entirely the nonequilibrium transport properties 
of the IHS gas. 

While the Boltzmann and Enskog equations for the 
EHS model are mathematically more tractable than 
those for the true IHS system, their solutions remain a 
formidable task. In the context of conventional fluids, 
those difficulties have stimulated the proposal of kinetic 
models, the prototype of which is the well-known BGK 
model kinetic equation 0, ITl| . An interesting practical 
application of the relationship IHS^EHS is the possibil- 
ity of extending those kinetic models to the case of gran- 
ular gases in a straightforward way. This has allowed us 
to recover simplified versions of models previously pro- 
posed for dilute and dense single gases In the case 
of inelastic mixtures, we are not aware of any previous 
proposal of a kinetic model, despite the fact that several 
kinetic models exist for elastic mixtures 01 • O ur method 
permits the construction of extended kinetic models for 
inelastic mixtures, as will be worked out elsewhere |39f . 
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factor aij, i.e., g'- a = — a^g ■ <r, and the kinetic energy 
decreases by a factor proportional to 1 — a? - , namely 

/2 ,2 2 2 / -^s2/i 2 \ miUlj 

rriiV +rrijVi —rriiV — m-iUi = — (g-er) (1— a,-.) . 

J J v ' v 1 ' mi + rrij 

(A4) 

From the velocity distribution function fi (v) of species 
i one can define the number density 

m = j dv/i(v), (A5) 

the mass density pi = m^rii, and the average velocity 

IK = - / dvv/ 4 (v) (A6) 
n% J 

of species i. The associated global quantities are the 
total number density n = Uj, the total mass density 
p = ^ZjPi, and the (barycentric) flow velocity 

u=-V / 5 1 u t . (A7) 

Pi 

The granular temperature T of the mixture is defined by 
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APPENDIX A: MIXTURES 



T = I J2 n ^ = ? / dv V2 f^), (A8) 



1. General properties 



where V = v — u is the peculiar velocity. In Eq. (|A8|) Tj is 
the partial granular temperature associated with species 
i. In general, equipartition of ener gy d oes not hold, even 
in homogeneous states, so Ti ^ T [20j. 

In the dilute limit, the distribution functions /i(v) 
obey a set of coupled Boltzmann equations, 

(9 t + v-V)/ l = ^J^ ) [/,,/ J ], (A9) 



In the case of a multi-component granular gas, the in- 
elastic collision between a sphere of species i (mass mi 
and diameter <Ji) and a sphere of species j (mass rrij and 
diameter <jj) is characterized by a coefficient of restitu- 
The direct and restituting collisions rules 



where the collision operator Jj^ [fi, fj] is given by 

J^ ) \fiJA = ^f 1 J dvi J d£e(g.5)( g .£) 
x h77*(v'%-K) - /iWiCvO] , (aio) 

where = (<Ji + <jj)/2. Every collision conserves the 
number of particles of each species, 



tion a 



Kl 



a, 



are given by 



v 

v', 



: v - pji (1 + atj) (g • 
vi + fiij (1 + a^) (g • (t)ct, 



= v 
vi 



(iji(l - 



Q-ij )(g • a-)a, 
-ar. 1 )(g-a)a, 



where 



fJ>ij 



(Al) 



(A2) 



(A3) 



Equations IjAip and (| A2|) are generalizations of Eqs. (|1.1|) 
and (|1.3|) . respectively. Again, the component of the 
post-collisional relative velocity along a is shrunk by a 



(AH) 



Moreover, the total momentum is conserved as well, i.e., 



J2m J dvvjg^LA,/,-] = o. (A12) 



However, the total energy is not conserved: 
JdvV 2 J^ ) [f i ,fj] = -±nT(, 



(A13) 
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what defines the collision rate £ of the mixture. 

In general, given an arbitrary function 4>(v), one can 
define its associated collision integral 



I 



4 rl / dv / dv i / d*e( g -*)(ga) 



x/i( v )/j( v i 



(A14) 



where in the last step we have performed a standard 
change of variables. Thus, Eqs. I|A11|) - (|A13I) can be 
rewritten as 

r(a»3-) 



I, 



'[i] = o, 



(A16) 



r(a«)r 



(A17) 



The collision integrals /^"'^[m^v] and /f" 1 ^ [m^y 2 ] 
cannot be evaluated exactly for arbitrary distribution 
functions fi and fj. This situation is analogous to the 
one taking place with Eq. I|2.8(l in the single gas case. 
Again, a reasonable estimate can be expected if the inte- 
grals are evaluated in the (multi-temperature) Gaussian 
approximation 

( m \ d/2 

/<(v) -> /i,o(v) = n t [—^\ exp {~m % V 2 /22}) , 

(A18) 

where we have restricted ourselves to the case = u 
in order to satisfy Eq. I|A6(I . When the app roximation 
l!8|) is inserted into Eq. (|A14|) one gets gfl E3 



TO,V 



0, 



(d- 



2)n,l 



m t (Tj 



'%]- 



OLi 



mjTi 



miTj 



(A19) 



(,A20) 



where 



4fi rf 



d-1 



2T\ 

mi 



1/2 



1 



mjTi 



3/2 



is an effective collision frequency of a particle of species 
i with particles of species j. In this Gaussian approxi- 
mation, the cooling rate defined by Eq. (|A17|) becomes 
C — * Co with 

Co = SI E (1 - 4) . (A22) 



4dnT 



where we have made use of the property piTvij = 
PjTji/ji. In the case of a single gas, Eqs. (|A21|) and (|A22|) 
reduce to Eqs. (|2.14|) and (|2.13|) . respectively. 



2. Model of frictional elastic hard spheres 



Once we have revised some of the basic properties of 
the collision operator J^" [fi, fj] for an IHS mixture, 
we are in conditions of proposing a minimal model for 
frictional EHS. In agreement with the philosophy behind 
Eq. H2.19JI . we write 



J, 



[fi i fj 



flij j\j ^ [fi i fj 



^•Kv-uO/,], 
(A23) 

where and Qj are to be determined by optimizing the 
agreement between the relevant properties of the true 



(A15) 

operator jj^^/j, fj] and those of the right-hand side of 
Eq. I|A23|) . For an arbitrary function </>(v), one has 

4T j) m - /v« m - 1 / dv (v - m) . ^/<(v). 

(A24) 

The above replacement satisfies Eq. (|A15(1 identically. In 
addition, it is consistent with Eq. (|A19|I in the Gaussian 
approximation (|A18|) . So far, j3ij and remain arbi- 
trary. Insertion of the Gaussian approximation (|A18|) 
into Eq. l|A24jl with 0(v) = niiV 2 yields 



^[miV 2 ] -> fl u (d + 2)n i T i i/ ij 



mjTi + niiTj 



(A25) 



where again we have restricted ourselves to the case Uj = 
u. Comparison between Eqs. (|A20|> and l|A25|) suggests 
the choices 



Ci: 



ft; 



d + 2 



(A26) 



(A27) 



Of course, other combinations of /3y and £y are in prin- 
ciple possible, but Eqs. I|A26|I and l|A27(l represent the 
simplest choice in the absence of mutual diffusion (i.e., 
with Ui = u). In the more general case Uj ^ u, the ex- 
pression for Cij is more complicated and will be reported 
elsewhere Q. 

Equation i|A25|) highlights that in general 
/( Qi j) [niiV 2 ] =/= because of two reasons. First, if 
species i and j have different mean kinetic energies (i.e., 
Tjj =/= Tj), then mutual collisions tend to "equilibrate" 
both partial temperatures. This equipartition effect, 
which is also present in the case of elastic collisions, 
is represented by the term f3ijI^[rriiV 2 ]. Moreover, 
even if T = Tj, one has I^[m t V 2 ] ^ due to the 
inelasticity of collisions, an effect that is accounted for 
by the term —QjolniT. While the former term can be 
either positive (T < Tj) or negative (T > Tj), the 
latter term is negative definite. Thus, Cij represents the 
cooling rate of species i due to collisions with particles 
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of species j. Since the relative decrease of energy after 
each i-j collision is proportional to 1 — a?^, it is quite 
natural that £y oc 1 — a?-. As in the single gas case, 
the parameter fyj measures the rate of i-j collisions of 
EHS relative to that of IHS. It is rather reinforcing that 
the choice (|A26(I is the natural extension of the choice 
(|3.31|) that we adopted for a single gas. While in the 
latter case we had to resort to the evaluation of the 
transport coefficients, the choice l|A26l) arises in the case 
of mixtures simply from the collisional energy transfer 
when Ti ^ Tj. 

3. Brownian limit 

Let us consider now the Brownian limit of a heavy im- 
purity particle (species 1) immersed in a bath of light par- 
ticles (species 2). In that case, the Boltzmann-Lorentz 
operator for inelastic collisions becomes the Fokker- 
Planck operator 

r(«i2)r f f i _ l+_ai2 ,(l)r, , , 1 7\ <9 2 

^12 [h,M - — ^ — J12 [h,m - 2^ 12 ^T^2'' 1 ' 

(A28) 

where 

4?lf„h]=^-(r + §-l) h . (A29) 
In the above equations, 

7U = ^^^^K/mO 1 / 2 {—\ 12 (A30) 
dy it \ mi J 

is the friction constant associated with elastic collisions 
and 

Ci2 = ^^(l-a? 2 ) (A31) 

is the cooling rate of the Brownian particle due to inelas- 
tic collisions with the bath particles. Equation (|A28|) is 
the exact Fokker-Plank limit of the inelastic Boltzmann- 
Lorentz operator. It is quite nice that Eq. (|A28(1 supports 
the choice (IA26|I as the collision rate factor. Moreover, 
Eq. (| A3 II for (12 agrees with the limit mi 3> m2 of the 
proposal l|A27jl . On the other hand, the exact Fokker- 
Planck equation (|A28|) differs from the one obtained from 
our model 1A23|) . More specifically, the latter results 
from the replacement 

C12 Ti d 2 (12 d 

2 mi ov z 2 ot 

While on the right-hand side of l|A32|) the role of mimick- 
ing the collisional cooling experienced by the Brownian 
particle is played by a deterministic frictional force, that 
role is played by a stochastic force on the left-hand side. 
On the other hand, both sides of l|A32|) agree in that they 
give vanishing contributions to the mass and momentum 



balance equations and yield the same contribution to the 
energy balance equation if Ui = u. In addition, they co- 
incide in the Gaussian approximation l|A18|) . Therefore, 
the model l|A23(l can be expected to behave reasonably 
well even in the extreme limit of a Brownian particle. 



4. Kinetic modeling 

As discussed in Sec. II VI the mapping 1A23(I allows one 
to transfer any given kinetic model 

jVifiJA^K^ (A33) 

for elastic mixtures [Tol ] into an equivalent model for in- 
elastic mixtures: 

jir\f, fj] - = p zjK y + & a . [(v - «*) /,] . 

(A34) 

Thus, one can construct in a straightforward way the 
inelastic versions of kinetic models such as Gross- 
Krook's Garzo-Santos-Brey's or Andries- 

Aoki-Perthame's ■ More details will be reported else- 
where (39) . 

APPENDIX B: DENSE GRANULAR GASES 

1. The Enskog equation for inelastic hard spheres 

So far, we have restricted ourselves to low density gran- 
ular gases for which the Boltzmann description seems 
appropriate. At higher densities the revised Enskog ki- 
netic theory, suitably generalized to inelastic collisions 
[hi|. provides a basis for analysis of granular flow. The 
Enskog kinetic equation reads 

(ft+v-V)/ = ^[/], (Bl) 

where 

4 a) [/]=cr d - 1 |dvi JdaQ(g-a)(g-a) 

x [a- a /a(r,r- < T J v",vl / )-/2(r > r + <T > v,vi)]. (B2) 

In this equation, er = era and 

/ 2 (r, r ± a, v, v x ) = *(r, r ± <x)/(r, v)/(r ± tr, Vi ) (B3) 

is the pre-collisional two-body distribution function in 
the molecular chaos approximation, x(r, r±er) being the 
equilibrium pair correlation function at contact as a func- 
tional of the density field. We use the notation [/] 

rather than jj^ [/, /] to remind that the nonlinear func- 
tional dependence of the Enskog collision operator on the 
velocity distribution function / is higher than bilinear, 
due to the presence of the correlation function \. 
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Taking velocity moments in (|B1|) one obtains again 
the hydrodynamic balance equations l|1.6|) - (|1.8|) . where 
the number density, flow velocity, and granular temper- 
ature are defined by Eqs. (|2.3I) - (|2.5[1 . respectively. The 
pressure tensor and the heat flux have both kinetic and 
collisional transfer contributions, namely P = Pk + P c , 
q = q*i + qc- The kinetic contributions Pk and are 
given by Eqs. (12. 9|) and (|2.1U|) . while the collisional trans- 
fer parts are [14( 



p(«) t/ ] 



i a) lf} = 



l±?-m<T d J dv J dvi J dae(g-a)(g-a) 2 

xerer / dA/ 2 (r - (1 - A)«t, r + Act, v, vi) , 

Jo 

(B4) 

J dv J dv x j d£6(g-S)(g-£) 2 



where now g = v /2T(r)/mC - ^2T(r + aj/md + 
u(r) — u(r + er). The local equilibrium cooling rate Co( r ) 
depends not only on the local values of the hydrodynamic 
fields at r but also on their values on a spherical surface 
of radius a around the point r. It also depends on the 
entire density field through x( r i r + <T )- If onc further 
neglects those dependencies, one gets the same result as 
in the dilute limit, Eqs. (|2.13fl and (|2.14fl . except that 
the collision frequency vq is multiplied by \. 



2. Model of frictional elastic hard spheres 

The natural extension to the Enskog equation of the 
model (|2~T§|) is 



-ma 



4 a3 [/]^/34 x) [/] + %l--(v/) 



2 <9v 



(Bll) 



x (G • a)a [ dX f 2 (r - (1 - X)a, r + Act, v, Vl ), Insertion of (BTTJ into Eqs. flU and |BS| yields 
Jo 

(B5) V- P { c a} ->/3V- pW, (B12) 



where G = (V + Vi)/2 = (v + v x )/2 - u(r). The su- 
perscript (a) on and qi^ has been introduced to 
emphasize that both quantities depend explicitly on the 
coefficient of restitution. They also depend implicitly on 
a through their functional dependence on the velocity 
distribution function /. With that convention, we can 
write 



V • q( Q) + P( Q) : Vu -» (3 V • + PW : Vu 



(Q) [/] = ± ^qi 1) [/]- (B6) 



The cooling rate is found to be [Til IT^ | 

C = (i-« 2 )"" 



ma 



AdnT 
x/ 2 (r,r + <r,v,vi). 



dv / dvi / dae(g-?)(g-?) 3 

(B7) 



The divergences of the collisional transfer parts are sim- 
ply related to moments of the collision operator > 



, (B13) 



where in Eq. (|B 1 3I> we have approximated C ^ Co- 
Comparison with the exact results (|B6|I implies that 
j3 = (1 + a)/2. Again, this reinforces the choice (|3.31() 
made in the dilute case. 



3. Kinetic modeling 

A few years ago, Dufty et al. 0, 03 proposed the 
following BGK-like kinetic model for the elastic Enskog 
equation: 



v • p! q) = -» 



dVV4 Q) [/], 



4 1} [/] 



(B8) 



-MO (/ - fo) - 
'mV 2 



fo 



dT 



- 1 



nT 
V • q« 



VV : PW 

Pi 1 ' : Vu 



V q 



(«) = _rn 



J dVV 2 4 a) [f]-P^ :Vu-^nCT. 

(B9) 

As in the dilute reasonable estimate of the cool- 

ing rate C(r, t) can be expected if the velocity distribution 
function is replaced by its local equilibrium approxima- 
tion H2.11fl . The resulting cooling rate is 



-AW : D(V) -BW • S(V) 



(B14) 



where 



A« = i J dvD(V)jW[/„], (B15) 



Co(r) = (I -a 2 ) 



! 



Ait d dT{v) 



dC / dC ie 



-c^-ct 



B (D 



2m 



(d + 2)T 2 



dvS(V)4 1} [/o] 



(B16) 



x / d5 ; 6(g-5 ; )(g-5 ; ) 3 x(r,r + cr)n(r + cr), 

(BIO) 



In Eqs. (Efll H BTg V D ( V ) and S ( V ) are S iven b y Eqs . 
(|3.9|) and i|3.11[l . respectively. Now, making use of I|B11|I . 
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the kinetic model (|B14|) is extended to inelastic collisions 
as 

1 



4%] 



^o(/-/o) + ||^-(V/) 



A 

nT 



VV : Pl Q) 



mV 



dT 



(V • q { c a) + P { c a} : Vu 



-A (q) : D(V) - B (q) • S(V) , (B17) 

where A( Q ) = ^(l + gjA' 1 ', B( Q ) = i(l + a)B«. The ki- 
netic model <|B17|1 for the inelastic Enskog equation turns 
out to be a simplified version of the model introduced in 
Ref. [3. 



2. Model of frictional elastic Maxwell particles 

In the case of frictional elastic Maxwell models 
(EMM), one makes the replacement l|2.19[l . The corre- 
sponding transport coefficients have the same forms as for 
EHS, so they are given again by Eqs. H3.17[) - I|3.19[l anc ^ 
(|3.3UI) . In the same spirit as before, if we want to optimize 
the agreement between the EMM and IMM transport co- 
efficients, the three possible choices for are 



_ r * c *_ (i + °0 2 

Pn ~ v f\ ^0 — ^ • 



-c 



(1 + a) 2 



(C7) 



(C8) 



APPENDIX C: INELASTIC MAXWELL MODELS 

1. Transport coefficients 

Now we return to a single dilute gas. In principle, 
the same philosophy behind Eq. (|2.19[1 can be applied 
to inelastic Maxwell models (IMM). In that case, the 
collision operator is [23, 

x/(r,v'/)-/(r,v)/(r,v 1 )]. (CI) 

This collision operator verifies again Eq. I|2.6|) but now 
the cooling rate is exactly given by ( — Coi Eq. i|2.13[l . 
The transport coefficients associated with the stress ten- 
sor and the heat flux have the structure of Eqs. (|3.3H 
(|3.5() . except that now the collision frequencies are |23| 



1 + a 



(1 + a)(d+l - a) 
2d 



±-l + h d + 8)(l-a) 



(C2) 



(C3) 



and the kurtosis of the homogeneous cooling state is 

6(1 -a) 2 



a hc = 

2 4d-7 + 3a(2-a)' 
Using the property 



J dvvjW[/ t ,/ hc ] = 



d + 2 



2d 



-fojt 



(C4) 



(C5) 



it is straightforward to prove that the self-diffusion coef- 
ficient is given by Eq. 1)3. 25|) . except that 



d+2 1 + a 
~Yd 2~' 



(C6) 



The Navier-Stokes transport coefficients characterizing 
an IMM mixture have been recently derived 0] . 



0d = 



2d 
d + 2 



D 



(C9) 



Interestingly, the shear viscosity and the thermal conduc- 
tivity routes give consistently the same expression for 
0, this expression being the square of the self-diffusion 
value. Since a\ c , which vanishes for EMM, does not ap- 
pear either in the shear viscosity or in the self-diffusion 
coefficient of IMM, the EMM model reproduces ex- 
actly those coefficients if = {3^ = (1 + a) 2 /4 and 
= fin = (1 + a)/2, respectively. However, even if 
= 0\ — (1 + a) 2 /4, the transport coefficients A and /i 
for IMM differ from those for EMM due to the explicitly 
appearance of a\ c in the former case. 

As Fig. shows, the choice = 0n = (1 + a)/2 de- 
scribes the behavior of 77, A, and fj, only at a rough qual- 
itative level. The alternative choice = 0„ = 0\ = 
(1 + a) 2 /4 strongly overestimates the diffusion coefficient 
and gives values for A and \i in reasonable agreement with 
the true IMM values only for small or moderate inelastic- 
ity. Therefore, the Navier-Stokes transport coefficients of 
IMM are very poorly described by the EMM model with 
a unique expression for 0, in contrast to what happens 
in the case IHS <-> EHS. This a consequence of the fact 
that the influence of inelasticity is much stronger in IMM 
than in IHS, as reflected on the kurtosis a^ and on the 
transport coefficients 77, A, /i, and D (compare Figs. |3| 
and[|JJ). As a matter of fact, the coefficients A and fj, di- 
verge at a = (4 — d)/3d for d — 2 and d = 3 23], thus 
indicating the failure of a hydrodynamic description for 
IMM with a < (4 - d) /3d. 

While, with the exception of the self-diffusion coeffi- 
cient |49| , the Navier-Stokes transport coefficients of IHS 
and IMM differ significantly [23|, those of frictional EHS 
and EMM agree each other [see Eqs. (|3~T7V(|3~lD|l and 
(|3.3U|) ]. provided the former are evaluated in the first So- 
nine approximation and the same value of is used in 
both approaches. As a matter of fact, hard spheres and 
Maxwell particles are known to exhibit similar rheologi- 
cal properties in the elastic case, even in states far from 
equilibrium 01 • This shows that the EMM system with 
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FIG. 6: (Color online) Plot of the (reduced) shear vis- 
cosity i)/(nT/i/|)), thermal conductivity \/(nT/mvo), trans- 
port coefficient [i/(T 2 /muo), and self-diffusion coefficient 
D/(T/mvo) for three-dimensional IMM ( — ) and the "equiv- 
alent" system of EMM with f3 = p v = P\ ( ), P = Pd (- • 

- • -), and P = 1 (-•■-••-) as functions of the coefficient of 
restitution a. Note that in the top panel the IHS curve and 
the EHS curve with p = p n = P\ are identical. The same 
happens in the bottom panel between the IHS curve and the 
EHS curve with p = p D . 



j3 = (1 + a)/ 2 can be usefully exploited as a model of 
IHS rather than as a model of IMM, namely 

1 + a _m , . Co d 
~2 



~-J 



EMM 



[/,/] 



(V/). (C10) 



Since the collision operator of EMM is mathematically 
more tractable than that of EHS, the mapping IHS — > 
EMM can be considered as intermediate between the 
mappings IHS — > EHS and IHS — > kinetic models. 



APPENDIX D: VELOCITY MOMENTS IN THE 
UNIFORM SHEAR FLOW 

In this Appendix we derive some expressions for the 
velocity moments predicted by the BGK and ES kinetic 
models for uniform shear flow. Let us focus on the 
isotropic moments 



{V 



2N\ 



dW^/(V). 



Insertion of the steady-state solution l|4.33[l yields 

r 2N \ _ /^oPr 



(V 



poo 

/ ds exp [- (/3z/ Pr + N( ) s] 
Jo 

x J dVf R (V) |A( S ) ■ V\ 2N , (D2) 



where Aij(s) = &y — as5i X 5j y . To carry out the velocity 
integral, we consider the diagonal representation of the 
matrix R _1 , namely 



S = ^ U .R- 1 .U-\ 
2 



(D3) 



where U is a unitary matrix whose expression will be 
omitted here. The eigenvalues are 



Si 



2Pr(/3 + Co*) 



2T 2(Pr/3 + Co*) ± (1 - Pr) [^/d(*(2f3 + dQ) T ^Co 

(D4) ' 



S, 



m Pr(/3 + C *) 
2T Pr/3 + Co ' 



(D5) 



where S3 is (d — 2)-fold degenerate. 

Using (|D3(1 . the velocity integral in Eq. I|D2|I becomes 



J dVf R (V) |A(s) • V| 2W = nn-^ 2 (dot S) 1/2 
x J dWe- S:WW |A(s) ■ IT 1 • W 



\2N 

where we have made the change of variables 

W = U • V. 

It is convenient to decompose the vector W as 
W = W x x + W y y + W ± , 

so that 



(D6) 

(D7) 
(D8) 



|A(s) • U" 1 • Wf = Wl + B xx (s)W 2 + 2B xy (s)W x W y 



where 



1 



Byy{s)Wl 

2~P + as^/dQ/2 



B Tr (s) = 1 + -a 2 s 2 — as 

2 VW+dQ 



1 



Byy(s) = 1 + —a s + as 



2P-a Sy /dQ/2 



Bxy(s) = as 



as^/pj2+ y/dQ 



(D9) 



(D10) 



(Dll) 



(D12) 



(Dl) Therefore, 
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p poo poo poo 

dVf R (V)\A( S )-V\ 2N = riir-^ 2 (detS) 1/2 tl d - 2 / dW x e-* w 2 dW y e-^ w y dW ± W^ 3 e~ s ^ 

J J — oo J— oq JO 



[Wl + B xx (s)W 2 + 2B xy (s)W x W y + B yy (s)W 2 ] 



N 



(D13) 



For any given value of N, the integrals over W x , W y , and Wj_ in Eq. (|D13jl . as well as the integral over s in Eq. (|D2J| . 
can be evaluated analytically. In particular, 



(V 4 



dT 2 



{ (2 + d)Pr 6 /3 8 + 2Pr 4 [l + 5(2 + d)Pr + (3 + d)Pr 2 ] ft 7 Q 



m2Pr/3(/3 + C *)'(Pr/3 + 2C *)° 
+Pr 2 [6d + 4(2 + d)Pr + 6(13 + 7d)¥i 2 + 20(3 + rf)Pr 3 + 3dPr 4 ] ft & Q 2 
+2Pr [3d + 4(1 + 4d)Pr + 2(34 + 25rf)Pr 2 + 40(3 + d)Pr 3 + 15dPr 4 ] ft 5 Q 3 
+2 [3d + 13dPr + 5(8 + 15d)Pr 2 + 80(3 + d)Pr 3 + 60dPr 4 ] f3 4 Q 4 
+2 [Ud - (16 - 41d)Pr + 2(122 + 39d)Pr 2 + 120dPr 3 ] (3 3 Q 5 

+2 [-16 + 28d + (112 + 15d)Pr + 120dPr 2 ] [3 2 Q 6 + 4 [8 + d(5 + 24Pr)] (3(* 7 + 18dQ 8 } , (D14) 

where use has been made of Eq. I|4.29|) . 

Non-isotropic velocity moments can be obtained in a similar way, although they are more complicated than the 
isotropic ones. In the case of the BGK model (Pr = 1), one has 



j3v Q 



ds exp [- {fivQ + (k + fc')Co) s] I dV (V x - asV y Y V^Vf / (V) 



(D15) 



After expanding (V x — asV y ) and carrying out the integrations over V and s, one finally gets 



(vX k - f vl k ') - {-if ( — 



2T 
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(25)! 
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